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Abstract. Gowdy waves, one of the standard 'apples with apples' tests, is 
proposed as a test-bed for constraint-preserving boundary conditions in the non- 
linear regime. As an illustration, energy-constraint preservation is separately 
tested in the Z4 framework. Both algebraic conditions, derived from energy 
estimates, and derivative conditions, deduced from the constraint-propagation 
system, are considered. The numerical errors at the boundary are of the same 
order than those at the interior points. 

Constraint-preserving boundary conditions is a very active research topic in 
Numerical Relativity [TJ [3J. During this decade, many conditions have been 
proposed, adapted in each case to some specific evolution formalism: Fritelli-Reula @j, 
Friedrich-Nagy 0, KST [fUE], Z4 [8], Generalized-Harmonic P UHl [H [2] , or BSSN [3j. 
Cross-comparison among different evolution formalisms has been carried out ('apples 
with apples' initiative [111 112]). But only periodic boundary conditions have been 
considered up to now. 

We endorse some recent claims (by Winicour and others) that the cross- 
comparison effort should be extended to the boundaries treatment. In this paper, 
we show that Gowdy waves [13], one of the 'apples with apples' tests, is suitable for 
boundary conditions cross-comparison in the non-linear regime. As an illustration, we 
test separately the energy-constraint preservation in the Z4 framework. We compare 
algebraic conditions, derived from energy estimates, with derivative conditions, 
deduced from the constraint-propagation system. The resulting numerical errors at 
the boundary are of the same order-of-magnitude than those at interior points. 

1. The Gowdy waves metric 

Let us consider the Gowdy solution |13j , which describes a space-time containing plane 
polarized gravitational waves. The line element can be written as 

ds 2 = t- 1 ' 2 e Q / 2 {-At 2 + dz 2 ) + t (e p dx 2 + e~ p dy 2 ) (1) 

where the quantities Q and P are functions of t and z only, and periodic in z. 
The initial slice t — to is usually chosen so that the simulations can start with an 
homogeneous lapse. 

Let us now perform the following time coordinate transformation 

* = t e- r ' T \ (2) 
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so that the expanding line element (|TJ) is seen in the new time coordinate t as 
collapsing towards the t = singularity, which is approached only in the limit r — > oo. 
This "singularity avoidance" property of the r coordinate is due to the fact that the 
resulting slicing by r = constant surfaces is harmonic |14j . We will run our simulations 
in normal coordinates, starting with a constant lapse ao = 1 at r = (t = to). 

Standard cross-comparison tests [111 I12j are currently done with periodic 
boundary conditions. But one gets basically the same results by setting up algebraic 
boundary conditions, which take advantage of the symmetries of the Gowdy line 
element. For a rectangular grid, planar symmetry allows trivial boundary conditions 
along the x and y directions. Also, allowing for the fact that the z dependence in (TTJ) 
is only through cos(27rz) , one can set reflecting boundary conditions for the interval 
< z < 1 . In this way, the Gowdy waves metric is obtained as a sort of stationary 
gravitational wave in a cavity with perfectly reflecting walls. 

We can then set up a full set of algebraic boundary conditions, which are 
consistent with the Gowdy line element Jl} for all times. This opens the door 
to a selective testing procedure, where one could for instance try some constraint- 
preserving condition for the longitudinal and transverse-trace modes, while keeping 
the exact condition for the transverse traceless ones. Or, as we will do below, testing 
just some energy-constraint preserving boundary conditions while dealing with all the 
remaining modes in an exact way. 



2. Characteristic decomposition 

We will consider here the the first-order version in normal coordinates, as described 
in refs. P~5J[II>]. For further convenience, we will recombine the basic first-order fields 
(Kij, Dijk, Aj, 0, Zi) in the following way: 

Hy = K ij -(trK-Q)>y ij Vi = j rs (D irs - D ris ) - Z h (3) 

IHjh = D ljk - {l rs D irs - Vi) j 3 k Wi=Ai- Y s D trs +2Vi (4) 

so that the new basis is (IF,-, fiijk, Wi, O, Vi). Note that the vector Z{ can be 
recovered easily from this new basis as 

Zi = -H k ik . (5) 

In order to compute the characteristic matrix, we will consider the standard form 
of (the principal part of) the evolution system as follows 

fit u + a a„F'(u) = ••• , (6) 

where u stands for the array of dynamical fields and F™ is the array of fluxes along 
the direction given by the unit vector n. With this choice of basic dynamical fields, 
the principal part of the evolution system gets a very simple form in the harmonic 
slicing case: 

F n {W l ) = F n (G) = V n F n (Vi) = m 6 (7) 
F n (IF, ) = X nij F n (fi kij ) = n k % (8) 

where the index n means a projection along m, and we have noted for short 

Kij = Mnij + n^Wj) - W n jij , (9) 

where round brackets denote index symmetrization. 
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We can now identify the constraint modes, by looking at the Fluxes of 9 and Zj 
in the array ([7][8]). It follows from that the energy-constraint modes are given by 
the pair 

e ± = e ± v n (io) 



with propagation speed ±a. Also, allowing for (|5I8[) . we can easily recover the flux of 
Zi: 

F n (Z l ) = -U n i (11) 

so that we can identify the momentum-constraint modes with the three pairs 

Mf = n m ± X nm , (12) 

with propagation speed ±a. Note that, allowing for ([3]), the longitudinal component 
n„„ does correspond with the transverse-trace component of the extrinsic curvature 
Kij. We give now the remaining modes: the fully transverse ones, with propagation 
speed ±a, 

Tab — Has ± KiAB (13) 

(the capital indices denote a projection orthogonal to n^), and the standing modes 
(zero propagation speed): 

Wi , V A , HAij ■ (14) 

Note that the standing modes (fTi|) . the energy modes (TIT)]) and the transverse 
momentum modes actually vanish for the Gowdy line element |T]). 

3. Energy-constraint preserving boundary conditions 

In refs. 15, 16J, the system above was shown to be symmetric hyperbolic, by providing 
a suitable energy estimate. We can rewrite it here as 

II,, 11" + X kl3 X kzj + 9 2 + V k V k + W k W k (15) 

This leads to the following sufficient condition for stability 

(IT^ \ nij + 9 V n ) |s > (16) 

where £ stands for the boundary surface (n being here the outward normal). 

Let us consider for instance the boundary at z — 1 . We can enforce there the 
partial set of exact (reflection) boundary conditions A n y — , so that the requirement 
(flljj) reduces to 

(6 V n ) |s > (17) 

which can be used for a separate test of energy-constraint preserving boundary 
conditions. 

As an illustration, we will test two such conditions. The first one is given in the 
form of a logical gate: 

(e v n ) |s < o => e |s = o (is) 

(9-gate), so that it only acts when condition (|17p is violated. The second one is an 
advection equation 

9 t 8|s =-ad n e-r)G, (19) 
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where we have included a suitable damping term. Note that the principal part of 
the constraint-preserving condition (|19p coincides with the 'maximal dissipation' one, 
dt E~ = , which is not constraint-preserving in the generic case. Condition (fT9|) 
can be understood as a sort of maximal dissipation condition for the energy-constraint 
evolution equation 



A 6 = 



(20) 



which follows from (the time component of) the covariant divergence of the Z4 field 
equations [T5]. 
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Figure 1. profiles along the z direction for a Gowdy waves simulation ending 
at r = 100. From bottom to top, results for the pure advection condition (119ft . 
damped advection (with rj = 0.2), 0-gate 1181 . and exact reflection (included here 
for comparison). 

We show in Fig. [1] our results for some numerical simulations ending at r = 100. 
We have included in the plot the exact (reflection) results for comparison. It is 
clear that the 0-gate condition gets very close to (actually slightly better than) the 
exact result in this case. The pure advection condition (|19p is off by half-an-order 
of magnitude. However, a suitable damping term (we have taken rj = 2) greatly 
improves this, leading in this case to even less error than the pure reflection condition. 
Note that we are showing here the Q profiles, giving the cumulated energy-constraint 
deviation. The average energy-constraint error in these strong-field simulations is 
actually smaller. 
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